/* Replication File for Doug and Gonzalo -Yearly */


/*Replication and Test File*/
/*READ BEFORE EXECUTION

1. please note that several packages are required to execute
the following analysis:

twowayfeweights
clusteff
boottest

2. The replication and additional tests rely on data from the Davis
and Wolfram replication file. Please download this file, execute 
their main-analysis.do file and then execute the following code
from the main-results folder, making sure to change the working
directory below to the correct location
*/

/*Stata Set-up*/
#delimit;
set matsize 10000;
set memory 2500m;
sysuse auto, clear;
set scheme s2color;
cd "\\babylon\geyerlab\jmaier\Research_Projects\Doug_Project\Steigerwald_VazquezBare_Maier2020\Davis2012\Data_Library_for_Davis_and_Wolfram_2012_AEJApp_2012_0052\main-results";



***TABLE 1: replication results with yearly balanced panel***
/* Create Balanced Yearly Panel */
#delimit;
use temp1.dta, clear;
drop if year < 1990;
by reactorid, sort: gen new_var = _N;
drop if new_var < 240;
drop if reactor_name=="Salem 1";
drop if reactor_name=="La Salle 1";
drop if reactor_name=="La Salle 2";
drop if reactor_name=="Clinton";
drop if reactor_name=="Millstone 2";
drop if reactor_name=="Millstone 3";
drop if reactor_name=="Crystal River 3";
drop if reactor_name=="Salem 2";
drop if reactor_name=="D.C. Cook 1";
drop if reactor_name=="D.C. Cook 2";
collapse(mean) power divested plantid age2 age3 capacity yeardivest, by(year reactorid);
replace divested = 1 if divested > 0;
save year_temp1.dta, replace;


/*Replication of Table 2 From Davis and Wolfram with a balanced yearly panel*/
#delimit;
use year_temp1.dta, clear;
quietly areg power divested, absorb(year) cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);
quietly areg power divested i.reactorid, absorb(year) cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);
quietly areg power divested i.reactorid age2 age3, absorb(year) cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);
quietly areg power divested i.reactorid age2 age3 [weight=capacity], absorb(year) cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);


/*Estimating Effective Number of Clusters*/
#delimit;
use year_temp1.dta, clear;
sum plantid;
sort plantid;
clusteff divested i.reactorid i.year, cluster(plantid) test(divested);


/*Wild boostrap standard error*/
#delimit;
use year_temp1.dta, clear;
regress power divested i.reactorid i.year, cluster(plantid);
boottest divested=0;



*** Figure 3: Newly Treated Reactors***

/* Numer of reactors treated in each period*/
#delimit;
use year_temp1.dta, clear;
gen DV = 0;
drop if divested == 0;
bysort reactorid (year): replace DV=1 if _n==1;
drop if DV ==0;
collapse(sum) divested, by(year);
rename divested new_divested;
save year_temp2.dta, replace;
use year_temp1.dta, clear;
collapse(sum) divested, by(year);
merge 1:1 year using year_temp2.dta;
replace new_divested = 0 if missing(new_divested);
save year_temp3.dta, replace;

/* graph new reactors treated */
#delimit;
use year_temp3.dta, clear;
graph bar new_divested, 
over(year, label(angle(45) labsize(vsmall)))
bar(1, color(gray))
ytitle("Number of Newly Divested Reactors")
graphregion(style(none) color(gs16))
ylabel(, nogrid);
graph export figures/newlytreatedtimeline_yearly.pdf, replace;


*** TABLE 2; WEIGHTS ***
/*weights*/
#delimit;
use year_temp3.dta, clear;
gen proportion_new = new_divested / 89;

/*Must manually calculate weights using equation and save as csv */

#delimit;
import delimited "weights_fromR_final.csv", clear;
drop v1;
gen id = _n;
reshape long group_, i(id) j(pid);
rename group_ weight_list;
drop if weight_list == 0;
histogram weight_list, frequency bin(50);
graph export figures/weight_hist.pdf, replace;
sum weight_list;
sum weight_list,d;



***For the next sets of analyzes we must identify the groups of interest (groups being those treated in the same year)***
/*Create Treatment Groups by Year*/
#delimit;
use year_temp1.dta, clear;
keep if divested > 0;
sort reactorid year;
bysort reactorid(year): keep if _n == 1;
egen group_id = group(year);
keep reactorid group_id;
save year_treatment_groups.dta, replace;

/*Merge data of groups*/
#delimit;
use year_temp1.dta, clear;
merge m:1 reactorid using year_treatment_groups;
replace group_id = 0 if group_id == .;
/* Save Data*/
save year_temp4.dta, replace;


/*Figure out number of treated groups in each period*/
#delimit;
use year_temp4.dta, clear;
gen DV = 0;
drop if divested == 0;
bysort reactorid (year): replace DV=1 if _n==1;
drop if DV ==0;
collapse(sum) group_id divested, by(year);


***Identified group 2 and group 3 as of interest do to large treatment group***

*** TABLE 3 ***

**GROUP 2**
/*Two Group Diff in Diff - G2*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==2); 
sum reactorid;
quietly reg power divested i.reactorid i.year, cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);

/*Estimating Effective Number of Clusters - G2*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==2); 
clusteff divested i.reactorid i.year, cluster(plantid) test(divested);

/*Wild boostrap standard error -G2*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==2) | (group_id==3);
xtreg power divested i.year if group_id!=3, fe i(reactorid) vce(cluster plantid);
boottest divested=0;

** GROUP 3 **
/*Two Group Diff in Diff -G3*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==3); 
sum reactorid;
quietly reg power divested i.reactorid i.year, cluster(plantid); lincom divested; display round(e(r2),.01); display e(N);

/*Estimating Effective Number of Clusters - G3*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==3); 
clusteff divested i.reactorid i.year, cluster(plantid) test(divested);


/*Wild boostrap standard error group - G3*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==2) | (group_id==3);
xtreg power divested i.year if group_id!=2, fe i(reactorid) vce(cluster plantid);
boottest divested=0;


***TABLE 4: Interacted Regression ***
#delimit;
use year_temp4.dta, clear;
gen control=0;
replace control=1 if (group_id==0);
gen entry2000=0;
replace entry2000 = 1 if (group_id==2) | (group_id==0);
gen entry2001=0;
replace entry2001 = 1 if (group_id ==3) | (group_id==0);
gen entry_2000_d = divested * entry2000;
keep if (group_id == 0)| (group_id==2) | (group_id ==3);
sum reactorid;
xtreg power entry_2000_d  divested i.year, fe i(reactorid) vce(cluster plantid);
boottest entry_2000_d  == 0;

***FIGURE 4: trends by group***
	  
/*Graph Power over Time by Group 2*/
#delimit;
use year_temp4.dta, clear;
collapse(mean) power, by(group year);
*drop if power < 0;
graph twoway (line power  year if group_id==0) ///
      (line power year if group_id==2), ///
	  legend(off)
	  ytitle("Average Power Output")
	  xline(2000, lcolor(black) lpattern(dash))
	  title("Trends by Treatment Group")
	  graphregion(style(none) color(gs16))
	  ylabel(, nogrid);
	  graph export figures/twogrouptrends_year2000.pdf, replace;
	  
	  
/*Graph Power over Time by Group 3*/
#delimit;
use year_temp4.dta, clear;
collapse(mean) power, by(group year);
*drop if power < 0;
graph twoway (line power  year if group_id==0) ///
      (line power year if group_id==3), ///
	  legend(off)
	  ytitle("Average Power Output")
	  xline(2001, lcolor(black) lpattern(dash))
	  title("Trends by Treatment Group")
	  graphregion(style(none) color(gs16))
	  ylabel(, nogrid);
	  graph export figures/twogrouptrends_year2001.pdf, replace;

	  
***FIGURE 5: dyanmic treament effects ***  

/*Two Group Event Study group 2*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==2); 
collapse(mean) divested yeardivest plantid power, by (year reactorid);
gen eventtime = year - yeardivest;
gen placebotime = -1 * eventtime;
replace eventtime = . if eventtime < 0;
replace placebotime = . if placebotime < 1;
/*Generate Treatment Dummies*/
levelsof eventtime, local(levels);
foreach l of local levels {;
gen B`l' = (eventtime==`l');
};
/*Generate Placebo Dummies*/
levelsof placebotime, local(levels);
foreach l of local levels {;
gen P`l' = (placebotime==`l');
};
/*Event Study Regression*/
#delimit;
reg power B0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P9 B1 B2 B3 B4 B5 B6 B7 i.reactorid i.year, vce(bootstrap, reps(1000) nodrop);
coefplot, keep(B0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P9 B1 B2 B3 B4 B5 B6 B7) vertical drop(_cons) yline(0, lcolor(black) lpattern(dot)) xline(10, lcolor(black) lpattern(dash)) graphregion(style(none) color(gs16))
	  ylabel(, nogrid)
xtitle("Event Time")
ytitle("Estimated Coefficient and Confidence Interval") 
title("Dynamic Effect of Divestiture on Reactor Performance")
order(P10 P9 P8 P7 P6 P5 P4 P3 P2 P1 B*);
graph export figures/twogroupeventstudy_yearly2.pdf, replace;
test P1 = P2 = P3 = P4 = P5 = P6 = P7 = P8 = P9 = 0;
est table, star(.05 .01 .001) keep(P* B*);


/*Two Group Event Study group 3*/
#delimit;
use year_temp4.dta, clear;
keep if (group_id == 0) | (group_id==3); 
collapse(mean) divested yeardivest plantid power, by (year reactorid);
gen eventtime = year - yeardivest;
gen placebotime = -1 * eventtime;
replace eventtime = . if eventtime < 0;
replace placebotime = . if placebotime < 1;
/*Generate Treatment Dummies*/
levelsof eventtime, local(levels);
foreach l of local levels {;
gen B`l' = (eventtime==`l');
};
/*Generate Placebo Dummies*/
levelsof placebotime, local(levels);
foreach l of local levels {;
gen P`l' = (placebotime==`l');
};
/*Event Study Regression*/
#delimit;
reg power B0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P9 B1 B2 B3 B4 B5 B6 B7 i.reactorid i.year,   vce(bootstrap, reps(1000) nodrop);
coefplot, keep(B0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P9 B1 B2 B3 B4 B5 B6 B7) vertical drop(_cons) yline(0, lcolor(black) lpattern(dot)) xline(10, lcolor(black) lpattern(dash)) graphregion(style(none) color(gs16))
	  ylabel(, nogrid)
xtitle("Event Time")
ytitle("Estimated Coefficient and Confidence Interval")
title("Dynamic Effect of Divestiture on Reactor Performance")
order(P10 P9 P8 P7 P6 P5 P4 P3 P2 P1 B*);
graph export figures/twogroupeventstudy_yearly3.pdf, replace;
test P1 = P2 = P3 = P4 = P5 = P6 = P7 = P8 = P9 = 0;
est table, star(.05 .01 .001) keep(P* B*);




